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Abstract 

A thermodynamically consistent damage model is proposed for the simulation of 
progressive delamination in composite materials under variable-mode ratio. The 
model is formulated in the context of Damage Mechanics. A novel constitutive 
equation is developed to model the initiation and propagation of delamination. A 
delamination initiation criterion is proposed to assure that the formulation can 
account for changes in the loading mode in a thermodynamically consistent way. 
The formulation accounts for crack closure effects to avoid interfacial penetration of 
two adjacent layers after complete decohesion. The model is implemented in a finite 
element formulation, and the numerical predictions are compared with experimental 
results obtained in both composite test specimens and structural components. 
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1 Introduction 


Structural collapse in a composite structure is caused by the evolution of 
different types of damage mechanisms, such as matrix transverse cracking, 
fibre fracture and delamination. The particular damage modes depend upon 
loading, geometry, lay-up and stacking sequence. 

Delamination is one of the most common types of damage in laminated fibre- 
reinforced composites due to their relatively weak interlaminar strengths. De- 
lamination may arise under various circumstances, such as in the case of trans- 
verse concentrated loads caused by low velocity impacts. This damage mode 
is particularly important for the structural integrity of composite structures 
because it is difficult to detect during inspection. Furthermore, delamination 
causes a drastic reduction of the bending stiffness of a composite structure 
and, when compressive loads are present, promotes local buckling. 

When other material non-linearities can be neglected, methods based on Lin- 
ear Elastic Fracture Mechanics (LEFM) have been proven to be effective in 
predicting delamination growth. However, LEFM cannot be applied without 
an initial crack. In some situations, methods combining a stress analysis with 
a characteristic distance have been applied to predict the initiation of de- 
lamination [l]-[2] . After delamination onset, LEFM can be used to predict 
delamination growth [3]- [4], Techniques such as virtual crack closure tech- 
nique (VCCT) [5]-[9], J-integral method [10], virtual crack extension [11] and 
stiffness derivative [12] have often been used to predict delamination growth. 
These techniques are used to calculate the components of the energy release 
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rate. Delamination growth is predicted when a combination of the components 
of the energy release rate is equal to, or greater than, a critical value [13]. 

However, difficulties are also encountered when these techniques are imple- 
mented using finite element codes. The calculation of fracture parameters, 
e.g. stress intensity factors or energy release rates, requires nodal variable and 
topological information from the nodes ahead and behind the crack front. Such 
calculations can be done with some effort for a stationary crack, but can be 
extremely difficult when progressive crack propagation is involved. 

Another approach to the numerical simulation of the delamination can be 
developed within the framework of Damage Mechanics. Models formulated 
using Damage Mechanics are based on the concept of the cohesive crack model: 
a cohesive damage zone or softening plasticity is developed near the crack 
front. The origin of the cohesive crack model goes back to Dugdale [14] who 
introduced the concept that stresses in the material are limited by the yield 
stress and that a thin plastic is generated in front of the notch. Barenblatt [15] 
introduced cohesive forces on a molecular scale in order to solve the problem 
of equilibrium in elastic bodies with cracks. Hillerborg et ah [16] proposed 
a model similar to the Barenblatt model, but where the concept of tensile 
strength was introduced. Hillerborg’s model allowed for existing cracks to grow 
and, even more importantly, also allowed for the initiation of new cracks. 

Cohesive damage zone models relate tractions to displacement jumps at an 
interface where a crack may occur. Damage initiation is related to the in- 
terfacial strength, i.e., the maximum traction on the traction-displacement 
jump relation. When the area under the traction-displacement jump relation 
is equal to the fracture toughness, the traction is reduced to zero and new crack 
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surfaces are formed. The advantages of cohesive zone models are their sim- 
plicity and the unification of crack initiation and growth within one model. 
Moreover, cohesive zone formulations can also be easily implemented in fi- 
nite element codes using decohesion elements [IT]- [25] . Although the cohesive 
damage models cannot be considered non-local damage models [26], they al- 
low a mesh-independent representation of material softening, provided that 
the mesh is sufficiently refined. 

In the formulation of the cohesive models, it is important to control the en- 
ergy dissipation during delamination growth in order to avoid the restoration 
of the cohesive state, i.e., it is necessary to assure that the model satisfies 
the Clausius- Duhem inequality. There are some models in the literature that 
can be used under constant mixed-mode conditions [21] , [23]- [24] , [30]- [35]. 
However, the models proposed generally do not satisfy the Clausius- Duhem 
inequality under variable-mode loading situation. Most of the models cited 
above define the damage threshold parameter as the maximum displacement. 
This assumption may lead to the violation of the Clausius-Duhem inequality 
when the crack grows in a varying mode. 

The restoration of the cohesive state is illustrated in Figure 1. This Figure 
represents the traction (r)-relative displacement (A) relation for two different 
mode ratios, Gu/ ( Gn + Gj) = A (Figure 1 a)) and Guf ( Gu + Gj) = B 
(Figure 1 b)), where Gj and Gjj are the components of the energy release 
rate. If the mode ratio changes from A to B during delamination growth, 
there is a restoration of the cohesive state. This effect is clearly inconsistent 
with the thermodynamical principles. 

[Figure 1 about here] 
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A damage model for the simulation of delamination under variable-mode is 
presented in this paper. A new delamination initiation criterion is proposed. 
The delamination onset criterion stems from the expression of the critical 
energy release rate for delamination propagation under mixed-mode loading 
proposed by Benzeggagh and Kenane [36]. The model is implemented in the 
implicit finite element code ABAQUS [54] by means of a user- written decohe- 
sion element. 

This paper is structured as follows: first, the formulation of the damage model 
for the simulation of delamination onset and growth model is presented. The 
finite element discretization of the boundary value is described. Finally, the nu- 
merical predictions are compared with experimental results obtained in com- 
posite test specimens and composite structural components. 


2 Model for delamination onset and propagation 

The boundary value problem, the kinematic equations, and the constitutive 
relations are presented for the formulation of the model for delamination onset 
and delamination propagation. 

2.1 Boundary value problem 

Consider a domain Q. as shown in Figure 2(a), containing a crack T c . The 
part of the crack on which a cohesive law is active is denoted by F f . 0 / ( and is 
called the fracture process zone (FPZ). 

[Figure 2 about here] 
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Prescribed tractions, t{, are imposed on the boundary P p, whereas prescribed 
displacements are imposed on P M . The stress held inside the domain, a t j , is 
related to the external loading and the closing tractions r+, rj in the cohesive 
zone through the equilibrium equations: 


a ijJ = 0 in hi 

(1) 

(jjjjij = ti on r p 

(2) 

rj = -rf = aijnj on Y coh 

(3) 


2.2 Kinematics of the interfacial surface 

To develop the necessary kinematic relationships, consider the crack T c shown 
in Figure 2(a) as part of a material discontinuity, T^, which divides the domain 
hi into two parts, 0+ and hi- (Figure 2(b)). 

The displacement jump across the material discontinuity [iq] , can be writ- 
ten as: 

M = u t ~ u 7 ( 4 ) 

where uf denotes the displacement of the points on the surface of the material 
discontinuity of the parts hl± of the domain. 

The fundamental problem introduced by the interfacial surface Y ( i is how to 
express the virtual displacement jumps associated to the surfaces r^± as a 
function of the virtual displacements. Consider a three-dimensional space with 
Cartesian coordinates Xj, i = 1,2,3 . Let the Cartesian coordinates xf de- 
scribe the position of the upper and lower surfaces F r/ + in the deformed con- 
figuration. Any material point on T d ± in the deformed configuration is related 
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to its undeformed configuration through: 


xf = Xi + uf 


( 5 ) 


where uf are the displacements with respect to the fixed Cartesian coordinate 
system. The coordinates x t of the midsurface can be written as [37]: 


Xi 




( 6 ) 


[Figure 3 about here] 


The components of the displacement jump vector are evaluated at the mid- 
surface r d , which is coincident with T f / in the undeformed configuration (see 
Figure 3). The midsurface coordinate gradients define the components of the 
two vectors, v v . and , that define the tangential plane at a given point, P: 



( 7 ) 

V U = Xi£ 

( 8 ) 


where r], £ are curvilinear coordinates on the surface F d . Although v v . and 
are generally not orthogonal to each other, their vector product defines a 
surface normal. Therefore, the local normal coordinate vector is obtained as: 


V,=V(X V, |[v £ X V, 


( 9 ) 


The tangential coordinates are then obtained as: 

v^vdW 1 

Vt = v n X V s 


( 10 ) 

( 11 ) 
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The components of v n , v s and v t represent the direction cosines of the local 
coordinate system in the global coordinate system at a material point P G T^. 
The director cosines define an orthogonal rotation tensor ; relating the local 
coordinate system to the fixed coordinate system. 

Using the rotation tensor, the normal and tangential components of the dis- 
placement jump tensor expressed in terms of the displacement held in global 
coordinates are: 

A m = 0 mi [tij] (12) 

where A m is the displacement jump tensor in the local coordinate system. 


2.3 Constitutive laws 

A constitutive law relating the cohesive tractions, Tj, to the displacement 
jump in the local coordinates, A*, is required for modeling the behavior of 
the material discontinuity. The constitutive laws in the material discontinuity 
may be formally written as: 

Tj = T (A i) (13) 

rj = D^A Z (14) 

where Z)U n is the constitutive tangent stiffness tensor. 

A new constitutive model relating the displacement jumps to the tractions, 
and based on Damage Mechanics is proposed. 

The delamination model proposed follows the general formulation of Contin- 
uum Damage Models proposed by Simo and Ju [38]- [39] and Mazars [40]. 
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The free energy per unit surface of the interface is defined as: 


V’ (A, d) = (1 — d) ■0° (A) (15) 

where d is a scalar damage variable, and ip° is a convex function in the dis- 
placement jump space defined as: 

V>° (A) = ^AiD^Aj * = 1,3; j = 1,3 (16) 


Negative values of A 3 do not have any physical meaning because interpen- 
etration is prevented by contact. Therefore, a modification of equation (15) 
is proposed to prevent interfacial penetration of the two adjacent layers after 
complete decohesion. The expression for the free energy proposed is: 

*/> (A, d) = (1 - d) */>° (Aj) - d*/>° ( 5 3i (-A 3 )) (17) 

where (•) is the MacAuley bracket defined as ( x ) = \ {x + |x|) and 5ij is the 
Kronecker delta. The constitutive equation for the interface is obtained by 
differentiating the free energy with respect to the displacement jumps: 

Ti = Wk = (1 _ d) ~ 6 D ‘i~ hi ( ” A3> (18) 


The undamaged stiffness tensor, Z)k, is defined as: 

Dl = SyK (19) 

where the scalar parameter K is a penalty stiffness. The constitutive equation 
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can be written in Voigt notation as: 


' 


f > 



Tl 


A, 


0 

T 2 

1 

Q_ 

A2 

> — d K < 

0 

r 3 


A3 

< > 


(-As) 

< > 


( 20 ) 


The energy dissipation during damage evolution, H , represented in Figure 4 
for single-mode loading, can be obtained from: 


= -^d>0 
dd ~ 


( 21 ) 


[Figure 4 about here] 


The model defined by equation (18) is fully determined if the value of the 
damage variable d is evaluated at every time step of the deformation process. 
For that purpose, it is necessary to define a suitable norm of the displace- 
ment jump tensor, a damage criterion, and a damage evolution law, as will be 
described in the following sections. 

2.3.1 Norm of the displacement jump tensor 

The norm of the displacement jump tensor is denoted as A and is also called 
equivalent displacement jump norm. It is used to compare different stages of 
the displacement jump state so that it is possible to define such concepts as 
‘loading’, ‘unloading’ and ‘reloading’. The equivalent displacement jump is a 
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non-negative and continuous function, defined as: 


A — sj (A3) 2 + (A^ ear ) 2 (22) 

where A3 is the displacement jump in mode I, i.e. , normal to midplane, and 
A s hear h the Euclidean norm of the displacement jump in mode II and in 
mode III: 

A shear = \j (Ai)" + (A2)" ( 23 ) 


2.3.2 Damage criterion 

The damage criterion is formulated in the displacement jump space. The form 
of this criterion is: 


F (A*, r*) := A* - r* < 0 Vt > 0 (24) 

where t indicates the actual time and r f is the damage threshold for the current 
time. If r° denotes the initial damage threshold, then F > r° at every point 
in time. Damage initiation is produced when the displacement jump norm, A, 
exceeds the initial damage threshold, r°, which is a material property. 

A fully equivalent expression for equation (24) that is more convenient for 
algorithmic treatment is [41]: 

F (A 4 , r*) := G (A*) - G (V) <0 Vt > 0 (25) 

where G(-) is a suitable monotonic scalar function ranging from 0 to 1. G(-) 
will define the evolution of the damage value, and will be presented in the 
following section. 
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2.3.3 Damage evolution law 


The evolution laws for the damage threshold and the damage variable must be 
defined in the damage model. These laws are defined by the rate expressions: 


r = /i 


(26) 


d 


. dF (A, r) 

d\ 


. dG (A) 


(27) 


where /i is a damage consistency parameter used to define loading-unloading 
conditions according to the Kuhn- Tucker relations: 


/} > 0 ; F(A t ,r*)<0 ; fiF (X^r 1 ) =0 (28) 


From the previous equations, it is easy to prove that the evolution of the 
internal variables can be integrated explicitly [38]: 


r f = max |r°, max A s | 0 < s < t 


(29) 


d* = 



(30) 


which fully describes evolution of the internal variables for any loading-unloading- 
reloading situation. The scalar function G (•) defines the evolution of the dam- 
age value. For a given mixed-mode ratio, f3, the function proposed here is 
defined as: 


G (A) 


A f (A - A 0 ) 
A (A / - A 0 ) 


(31) 


Equation (31) defines the damage evolution law by means of a bilinear con- 
stitutive equation (see Figure 5), where A 0 is the onset displacement jump, 
and it is equal to the initial damage threshold r°. The initial damage thresh- 
old is obtained from the formulation of the initial damage surface or initial 


12 



damage criterion. is the final displacement jump, and it is obtained from 
the formulation of the propagation surface or propagation criterion. 


[Figure 5 about here] 


It is therefore necessary to establish the delamination onset and propagation 
surfaces for the complete definition of the damage model. Delamination on- 
set and propagation surfaces and the damage evolution law fully define the 
constitutive equations. 

The constitutive equations for the interfacial surface are normally developed in 
a phenomenological way, i.e., satisfying empirical relations that are obtained 
using experimental results. There are several types of constitutive equations 
used in decohesion elements: Tvergaard and Hutchinson [42] proposed a trape- 
zoidal law, Cui an Wisnom [43] a perfectly plastic rule, Needleman first pro- 
posed a polynomial law, [28], and later an exponential law [29]. Goyal et ah 
[44] adopted Needleman’s exponential law to account for load reversal without 
restoration of the cohesive state. 

The law proposed here is a bilinear relation between the tractions and the 
displacement jumps [21], [24], [45]. The bilinear law is the most commonly used 
cohesive law due to its simplicity. One drawback of the bilinear law is that 
the traction-displacement jump relation is discontinuous at the peak value 
of the traction. The discontinuity in the traction-displacement jump relation 
can be avoided using continuous functions. However, even for such continuous 
functions, the discontinuity is unavoidable when modeling loading-unloading 
cycles. 
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For a given mixed-mode ratio, /?, defined as: 


P 


A 


shear 


^ shear (^ 3 ) 


(32) 


the bilinear constitutive equation is defined by a penalty parameter, K, the 
damage value, d, the mixed-mode damage initiation, A 0 , and the total decohe- 
sion parameter, . These last two values are given by the formulation of the 
onset and the propagation criterion which takes into account the interaction 
between different modes, and their value depends on the mixed-mode ratio (3. 
The penalty parameter K assures a stiff connection between two neighboring 
layers before delamination initiation. The penalty parameter should be large 
enough to provide a reasonable stiffness but small enough to avoid numeri- 
cal problems, such as spurious tractions oscillations [46], in a finite element 
analysis. 


Propagation criterion 

The criteria used to predict delamination propagation under mixed-mode load- 
ing conditions are usually established in terms of the components of the energy 
release rate and fracture toughness. It is assumed that when the energy re- 
lease rate, G, exceeds the critical value, the critical energy release rate G c , 
delamination grows. 

The most widely used criterion to predict delamination propagation under 
mixed-mode loading, the ” power law criterion” is normally established in 
terms of a linear or quadratic interaction between the energy release rates 
[48]. However, Camanho et al. [24] have shown that the expression proposed 
by Benzeggagh and Kenane [36] for the critical energy release rate for a mixed- 
mode ratio is more accurate for epoxy and PEEK composites. The expression 
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proposed by Benzeggagh and Kenane for the critical energy release rate G c is: 


G c = G Ic + (G IIc - G Ic ) (33) 

Delamination growth is produced when the total energy release rate G is 
greater or equal than the critical value G c : 

G > G c (34) 


The energy release rate under mixed-mode loading is G = Gi + G s h ea r where 
G shear = Gn + G[[[ is the energy release rate for shear loading proposed by 
Li [49], [50], 

The propagation surface in the displacement jump space is defined through the 
final displacements, which are obtained from the pure mode fracture toughness 
( Gic, Gnc, Giuc ) considering that the area under the traction-displacement 
jump curves is equal to the corresponding fracture toughness, i.e. : 

G c = ^KA°A f (35) 


Using equation (35) in equation (33) the propagation criterion is obtained in 
the displacement jump space as: 


A f = 3 3 


A° W A 

shear A§A l 


G„ 


Gt 


A 0 


(36) 


where Ag and A° shear are the pure mode onset displacement jumps and A{ 
and A{ hear are the pure mode final displacement jumps. It is necessary to 
obtain the ratio G Q^ ar to fully define the final displacement jump. For a given 


mixed-mode ratio, /?, the energy release rates are obtained from: 


G, = \k (A» (0) A{ (0) - A 3 A 3 <fi)) 


(37) 
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G shear ~ shear (P) A shear (P) ^ shear ^{hear (/O) 


( 38 ) 


where A° shear (/5) and (/?) are respectively the shear and normal displace- 
ment jump corresponding to the onset of softening under mixed- mode loading, 
A {hear (P) an d A 3 (/?) are the shear and normal displacement jump correspond- 
ing to the total decohesion under mixed-mode loading, and A s h ea r and A 3 are 
the components of the current displacement jump. 


From (32): 



W) = ^ (« Ay 


(39) 


ALar (, 3) = A l tfi) Ay 


(40) 


A -A 13 

aa shear aA3 ^ ^ 


(41) 

Using equations (39), 

(40), and (41) in (37) and (38), the ratio between 

G shear 
Gt 

can be established in 

terms of /?. Since the ratio G ^y ar is 

only a function of 

the mixed-mode ratio 

/?, henceforward this ratio is named 

as B: 



t -> G shear ft 


(42) 


G t 1 + 2 (3 2 - 2/1 



Initial damage surface 

Under pure mode I, mode II or mode III loading, delamination onset occurs 
when the corresponding interlaminar traction exceeds its respective maximum 
interfacial strength, Under mixed-mode loading, an interaction be- 

tween modes must be taken into account. Few models take into account the 
interaction of the traction components in the prediction of damage onset. The 
models that account for the interaction of the traction components are usually 
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based on Ye’s criterion [51], using a quadratic interaction of the tractions: 



However, experimental data for the initiation of delamination under mixed- 
mode is not readily available and, consequently, failure criteria that can predict 
the initiation have not been fully validated. 

The criterion for propagation is often formulated independently of the crite- 
rion for initiation. In this paper, a link between propagation and initiation 
is proposed. Since delamination is a fracture process, the initiation criterion 
proposed in this paper evolves from the propagation criterion and the damage 
evolution law. The isodamage surface for a damage value equal to 1 corre- 
sponds to the propagation surface obtained from equation (33). Then, the 
isodamage surface for a damage value equal to 0 is the initial damage surface. 
With these assumptions, the criterion for delamination initiation proposed 
here is: 

(r») 2 = (r 3 ) 2 + (ri) 2 + (r 2 ) 2 = M) 2 + ((t^J 2 - (r») 2 ) (44) 

In the displacement jump space, the criterion becomes: 

(A») 2 = (A 3 ) 2 + (AO 2 + (AO 2 = (A») 2 + ((AL„ r ) 2 - (a §) 2 ) B n (45) 

The initiation criterion developed here and summarized by equation (44) is 
compared with Ye’s criterion and with a maximum traction criterion that does 
not take into account interaction between the tractions. The surfaces obtained 
by the different criterion are represented in Figure 6. The values predicted by 
the new criterion are very close to Ye’s criterion, that has been successfully 
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used in previous investigations [24], 


[Figure 6 about here] 

The formulation presented assures a smooth transition for all mixed-mode ra- 
tios between the initial damage surface to the propagation surface through 
damage evolution. If the loading mode changes, the formulation presented 
avoids the restoration of the cohesive state and assures that the energy dissi- 
pation is always positive. 

The evolution of the damage surface from the damage initiation surface to the 
propagation surface is represented in Figure 7, for positive values of displace- 
ment jumps. 


[Figure 7 about here] 


2-4 Formulation of the constitutive tangent tensor 


The constitutive tangent tensor needs to be defined for the numerical imple- 
mentation of the proposed model. The constitutive tangent tensor is obtained 
from the differentiation of the secant equation (18): 


T i — DijAj SijK 


1 + S 3 j- 


A'>" 


Ajd 


(46) 


where is defined as: 


Dij hjj F 


1 - d 1 + 5, 


-A,)' 

3J A,- 


(47) 
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The evolution of the damage variable d only occurs for loading situations. 
Then, the evolution of the damage variable can be written as: 


d 


G(A) = ^lA ,r<A<A f 

0 , r > A or A-f < A 


(48) 


where the variation of the function G is obtained assuming that the variation 
of the final displacement jump and the onset displacement jump A 0 with 
the mixed-mode ratio [3 are not significant for the time increment taken: 


dG (A) _ A^A° 1 
~dX~ ~ A/ - A° A^ 


(49) 


The evolution of the displacement norm is obtained from equation (22): 


A 


d\ : _ A fc 
dA k k A 



(-A fc ) \ 

A fc J 


Afc 


(50) 


Using equations (48) through (50), equation (46) can be written as: 


Ti 


D tan A ■ 


(51) 




tan 


{Dij - K [f + 6. 


izM" 

3 J A, 


l + Ssi^lHAAj} ,r < X < A f 


= < 




, r > A or A-f < A 
(52) 


where if is a scalar value given by: 


A f A 0 1 
A f - A 0 A 1 


( 53 ) 
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3 Finite element discretization - computational model 

To transform the strong form of the boundary value problem into a weak form 

better suited for finite element computations, the velocities V{ must belong 

to the set U of the kinematically admissible velocity held which allows for 
discontinuous velocities across the boundary T^ of the delamination. 

The spaces for the test functions and trial functions are defined as: 

5vi(X)eUo, Uo = [5vi\Svi e C° (X) ,5 Vi = 0 on T,.} (54) 

Vi(X,t)eU , U = [vi\vi e C° (X) ,Vi = Vi on T„.} (55) 

The space of velocities in U are the kinematically admissible velocities or 
compatible velocities. The space U satisfies the continuity conditions required 
for compatibility and the displacement boundary conditions. 

Considering Figure 2, the weak form of the momentum equation is obtained 
as: 


£ 


' do; 


v 


+ phi SvidQ± = 0 VTj e U 


k 2 ± \ dxj j 

where bi are the body forces and p is the density of the material 


(56) 


Using the decomposition of the velocity gradient and the traction continuity 
condition, the weak form of the momentum equation in an updated Lagrangian 
formulation is obtained as: 



US [nj] clT d + J2 / hDijOijd f2± = 
n± 

VjQji ) dT±+Y i / 5vipbidkl± Mvi&U 

o J 


(57) 
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where U is the traction tensor, (Ji 3 is the Cauchy stress tensor, and Dij is the 
rate of deformation tensor. From equation (57), it is clear that the tractions 
occurring at the cohesive interface are work-conjugate with the displacement 
jumps. 

The discretization of the domain is performed by the discretization of the 
whole domain hi with standard volume elements. However, the surfaces sur- 
rounding the potential delamination F f / are discretized with decohesion ele- 
ments [24] . The discretized formulation is divided in the two domains consid- 
ering no formal coupling between the continuous and the discontinuous parts 
of the deformation in the expression for the free energy of the interface [53]. 

3. 1 Discretization of the interfacial surface 

The displacements and displacement gradients for the decohesion elements are 
approximated as: 

«ik = Kil (58) 

M la = KC ( 59 > 

with: 

N e K Ke r+ 

N e K = (60) 

-N e K KeT d 

y. 

where q e Ki is the displacement in the i direction of the K node of the element 
e, N ( k are standard Lagrangian shape functions [52], N K are Lagrangian shape 
functions defined for the decohesion elements [24], 

According to equation (58), the displacement held, Ui, and the undeformed 
material coordinate, X i} associated with the surfaces T d ± are interpolated as 
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follows: 


uf = N K qf a 

(61) 

X f = N K p ki 

(62) 


where are the nodal displacement vectors and pf a are the undeformed 
material nodal coordinate vector. Note that the values of p^ Ki and p^ can be 
different in the case that an initial crack exists. Using these equations, the 
material coordinates of the interfacial midsurface are: 

Xi = - N Ki (pfci + p Kl + + q Ki ) (63) 

The components of the two vectors that define the tangential plane can be 
written as: 

% = = N KijV - (p^ + p~ Ki + ql H + fe) (64) 

= U,? = N K;,PT 1 ( PKi + V~Ki + qki + Qm) ( 65 ) 

Using (59) and (12), the displacement jump can then be obtained in local 
coordinates as: 

= Q i m N f( q j = BimKqKi (66) 

The contribution of a decohesion element for the internal load vector is given 
by: 

f$ = [ r n B inK dT d (67) 

J r d 

with B mK = Q in N 

The softening nature of the decohesion element constitutive equation causes 
difficulties in obtaining a converged solution for the non-linear problem when 
using Newton-Raphson iterative method. In particular, quadratic convergence 
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is not assured because the residual vector is not continuously differentiable 
with respect to the nodal displacements. 

The tangent stiffness matrix stems from the linearization of the internal force 
vector and it is obtained using Taylor’s series expansion about the approxima- 
tion qxi [25] . Taking into account that the calculation of the geometric terms 
of the tangent stiffness matrix is computationally very intensive, these terms 
are neglected. The tangent stiffness matrix, K r ziK , for the decohesion element 
is therefore approximated as: 

KrZiK ~ [ BjrzDlfBinKdTd (68) 

JTd 

where is the material tangent stiffness matrix, or constitutive tangent 
tensor defined in 2.4. 


4 Comparison with experimental studies 

The formulation proposed here was implemented in the ABAQUS Finite Ele- 
ment code [54] as a user- written element subroutine (UEL). 

To verify the element under different loading conditions, the numerical pre- 
dictions were compared with experimental data obtained for composite test 
specimens and aircraft subcomponents. The double cantilever beam (DCB) 
test, the end notched flexure (ENF) test, and mixed-mode bending (MMB) 
tests in PEEK/AS4, a thermoplastic matrix composite material, were simu- 
lated. 

The debonding of a composite co-cured skin- stiffener subcomponent loaded 
under tension was simulated, and the numerical results were compared with 
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experimental data. 


4-1 Mode I, mode II and mixed-mode I and II delamination growth for a 
PEEK composite 

The most widely used specimen for mixed-mode fracture is the mixed-mode 
bending (MMB) specimen shown in Figure 8, which was proposed by Reeder 
and Crews [55] -[57]. 


[Figure 8 about here] 

The main advantages of the MMB test method are the possibility of using vir- 
tually the same specimen configuration as for mode I tests, and the capability 
of obtaining different mixed-mode ratios, ranging from pure mode I to pure 
mode II, by changing the length c of the loading lever shown in Figure 8. 

An 8-node decohesion element is used to simulate DCB, ENF and MMB tests 
in unidirectional AS4/PEEK carbon-fiber reinforced composite. The spec- 
imens simulated are 102-mm-long, 25.4-mm-wide, with two 1.56-mm-thick 
arms. The material properties are shown in Table 1, and a stiffness K = 10 6 
N/mrn 3 is used. 


[Table 1 about here] 

The experiments used to assess the accuracy of the model proposed were per- 
formed by Reeder [55]- [57]. The experimental tests were performed at different 
jP ratios, ranging from pure mode I loading to pure mode II loading. The 
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initial delamination length of the specimens (oo) and the mixed-mode fracture 
toughness obtained experimentally are shown in Table 2. 


[Table 2 about here] 


Models using 150 decohesion elements along the length of the specimens, and 
4 decohesion elements along the width, were created to simulate the ENF 
and MMB test cases. The initial size of the delamination is simulated by 
placing open decohesion elements along the length corresponding to the initial 
delamination of each specimen (see Table 2). These elements are capable of 
dealing with the contact conditions occurring for mode II or mixed-mode I 
and II loading, therefore avoiding interpenetration of the delamination faces. 
The model of the DCB test specimen uses 102 decohesion elements along the 
length of the specimen. 

The different Gn/Gi ratios are simulated by applying different loads at the 
middle and at the end of the test specimen. The analytical determination of 
the middle and end loads for each mode ratio is presented in [24], 

The experimental results relate the load to the displacement of the point 
of application of the load P in the lever (load-point displacement, Figure 
8). Since the lever is not simulated, it is necessary to determine the load- 
point displacement from the displacement at the end and at the middle of the 
specimen using the analytical procedure described in [24], 

The B-K parameter, rj = 2.284, is calculated by applying a least-squares fit to 
the experimental data shown in Table 2. 
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Figures 9 to 13 show the numerical predictions and the experimental data for 
all the test cases simulated. 


[Figure 9 about here] 


[Figure 10 about here] 


[Figure 11 about here] 


[Figure 12 about here] 


[Figure 13 about here] 


Table 3 shows the comparison between the predicted and experimentally de- 
termined maximum loads. 


[Table 3 about here] 


It can be concluded that a good agreement between the numerical predictions 
and the experimental results is obtained. The largest difference (—8.1%) cor- 
responds to the case of an MMB test specimen with = 20%. This fact 
is not surprising, since the largest difference between the fracture toughness 
experimentally measured and the one predicted using the B-K criterion occurs 
for = 20%. 
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4-2 Skin- stiffener co-cured structure 


Most composite components in aerospace structures are made of panels with 
co-cured or adhesively bonded frames and stiffeners. Testing of stiffened panels 
has shown that bond failure at the tip of the stiffener flange is a common failure 
mode. Comparatively simple specimens consisting of a stringer flange bonded 
onto a skin have been developed by Krueger et al. to study skin /stiffener 
debonding [58]. The configuration of the specimens studied by Krueger is 
shown in Figure 14. 


[Figure 14 about here] 

The specimens are 203 mm-long, 25.4 mm-wide. Both skin and flange were 
made from IM6/3501-6 graphite/epoxy prepreg tape with a nominal ply thick- 
ness of 0.188 mm. The skin lay-up consisting of 14 plies was (0°/45°/90°/- 
45 o /45°/-45 o /0°)s and the flange lay-up consisting of 10 plies was (45°/90% 
45°/0°/90°) s . 

The properties of the unidirectional graphite/epoxy and the properties of the 
interface reported in reference [58] are shown in Tables 4 and 5, respectively. 

[Table 4 about here] 

[Table 5 about here] 

The parameter for the B-K criterion is taken from test data for AS4/3501-625 
as 77=1.45, and a stiffness K = 10 6 N/nnn 3 is used. 
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To keep the modeling difficulties low and the approach applicable to larger 
problems, the model that was developed uses only two brick elements through 
the thickness of the skin, and another two through the flange. The complete 
model consists of 1,002 three-dimensional 8-node elements and 15,212 degrees 
of freedom. To prevent delamination from occurring at both ends of the flange 
simultaneously, model symmetry was reduced by modeling the tapered end of 
the flange with a refined mesh on one side and a coarser mesh on the other. 
Unlike the previous examples, this model does not contain any pre-existing 
delaminations. 

Residual thermal effects in the composite plies are simulated by performing 
a thermal analysis step before the mechanical loads are applied. The same 
coefficients of thermal expansion (an=-2.4xlCU 8 /°C and a22=3.7xl0 -5 /°C) 
are applied to the skin and the flange, and the temperature difference between 
the room and curing temperatures is AT=-157 °C. The flange has more 90° 
plies than 0° plies, and the skin is quasi-orthotropic, so it is expected that 
residual thermal stresses are present at their interface at room temperature. 

Deformed plots of the finite element model immediately before and after flange 
separation are shown in Figure 15. 

[Figure 15 about here] 


It can be observed that only the refined end of the flange separates. It is 
worth noticing that the debond growth is not symmetric across the width: the 
debond initiates on the left corner of the flange shown on the bottom left of 
Figure 15 due to the lack of symmetry introduced by the terminated plies at 
the flange tapered ends. This behavior was also observed in the experiments 
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[ 58 ], 


Figure 16 shows the load-extensometer measurement relation obtained in 4 
experiments and the corresponding numerical prediction. Debond is detected 
in the experiments by the discontinuities in the load-displacement relation. 

Table 6 compares the average of the measured debond loads with the numerical 
predictions. 


[Figure 16 about here] 

[Table 6 about here] 

It can be observed that good accuracy in the prediction of the debond loads 
is obtained. The predicted stiffness of the specimen is also in good agreement 
with the experimental data. The stiffening effect detected in the experiments, 
Figure 16, is due to the extensometer rotation as a result of specimen bending. 
Although specimen bending is properly represented by the numerical model, 
the extensometer measurement calculated from the numerical model does not 
account for the rotation of the extensometer. 


5 Concluding remarks 

A thermodynamically consistent model for the simulation of progressive de- 
lamination based on Damage Mechanics was presented. A constitutive equa- 
tion for the interface was derived from the free energy of the interface. The 
resulting damage model simulates delamination onset and delamination prop- 
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agation. The constitutive equation proposed uses a single scalar variable to 
track the damage at the interface under general loading conditions. A new 
initiation criterion that evolves from the Benzeggagh-Kenane propagation cri- 
terion has been developed to assure that the model accounts for changes in the 
loading mode in a thermodynamically consistent way and avoids restoration 
of the cohesive state. The damage model was implemented in a finite element 
model. The material properties required to define the element constitutive 
equations are the interlaminar fracture toughnesses, the penalty stiffness, and 
the strengths of the interface. In addition, a material parameter r/, which is 
determined from standard delamination tests, is required for the Benzeggagh- 
Kenane mode interaction law. 

Two examples were presented that test the accuracy of the method. In the first 
example, the simulations of the DCB, ENF and MMB tests represent cases 
of single-mode and mixed-model delamination. A composite skin-stiffener co- 
cured sub-component was also simulated, and the model predictions were com- 
pared with available experimental data. 

The examples analyzed are in good agreement with the test results, and they 
indicate that the proposed formulation can predict the strength of composite 
structures that exhibit progressive delamination. 
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Figures 


Loading mode A Loading mode B 

G,,/ (Gn+G|)=A G II /(G II +G,)=B 



Fig. 1. Restoration of the cohesive state for delanrination propagation under vari- 
able-mode ratio. 
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Fig. 2. Body fl crossed by a material discontinuity in the undeformed configu- 
ration. 



Undeformed configuration Deformed configuration 


Fig. 3. Interfacial surface deformation. 
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Fig. 4. Energy dissipation during damage evolution. 


Onset 



Fig. 5. A bilinear constitutive equation for the decohesion element for a mixed mode 
loading situation. 
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Fig. 6. Comparison between Ye’s criterion, a maximum traction criterion and the 
new proposed criterion. 



A 3 (mm) 


Fig. 7. Damage evolution surface in the relative displacements space. 
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Fig. 8. MMB test specimen. 



Fig. 9. Numerical and experimental results- pure mode I loading. 
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Fig. 10. Numerical and experimental results- mixed mode I and II loading with 
G///Gt=20%. 
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Fig. 11. Numerical and experimental results- mixed mode I and II loading with 
G h /G t =50%. 
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Fig. 12. Numerical and experimental results- mixed mode I and II loading with 
G///Gt=80%. 



Fig. 13. Numerical and experimental results- pure mode II loading. 
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Fig. 14. Skin-stiffener test specimen. 



Fig. 15. Skin-Stiffener debonding. 
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Fig. 16. Experimental and numerical load-extensometer displacement relations. 
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Table 1 


Tables 


Ply properties. 


Ell E'22 — E33 G12 — G13 G23 v 12 — ^13 


122.7 GPa 10.1 GPa 5.5 GPa 3.7 GPa 0.25 


^23 

Gic 

Guc 


r° 

r 3 

II 

OcN 

0.45 

0.969 kJ/rn 2 

1.719 kJ/rn 2 

80 MPa 100 MPa 

Table 2 

Experimental data. 

Gill Gt 

0% (DCB) 

20% 

50% 

80% 

100% (ENF) 

G c [kJ/rn 2 ] 

0.969 

1.103 

1.131 

1.376 

1.719 

ao [mm] 

32.9 

33.7 

34.1 

31.4 

39.2 
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Table 3 


Maximum loads. 


Gu/Gt 

Predicted [N] 

Experimental [N] 

Error (%) 

0% (DCB) 

152.4 

147.5 

3.4 

20% 

99.3 

108.1 

-8.1 

50% 

263.9 

275.3 

-4.2 

80% 

496.9 

518.7 

-4.2 

100% (ENF) 

697.1 

748.4 

-6.9 


Table 4 


Material properties for IM6-3501-6 unidirectional graphite epoxy. 


Ei (GPa) 

E 2 =E 3 (GPa) 

M 2 — PL3 ^23 Gl 2 — Gi 3 (GPa) 

G 23 (GPa) 

144.7 

9.6 

0.3 0.45 5.2 

3.4 
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Table 5 


Interface properties. 


G ic (Nnrnr x ) 

G uc (Nnrnr x ) 

N (MPa) 

S (MPa) r] 

0.075 

0.547 

61 

68 1.45 


Table 6 

Comparison between experimental and numerical results. 


Experimental (kN) 

Predicted (kN) 

Error (%) 

Flange debond load 22.7 

21.0 

-7.5 
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